Fatigue-life distribution (fatiguelife) — Birnbaum–Saunders model for fatigue failures#
The fatigue-life distribution (also known as the Birnbaum–Saunders distribution) is a positive, right-skewed continuous distribution used to model time-to-failure under cyclic fatigue.
A key characterization is a normalizing transform to a standard normal:
equivalently, a constructive representation:
Learning goals#
understand what
fatiguelifemodels and when it is usedwrite down the PDF/CDF and connect them to a standard normal
compute moments (mean/variance/skewness/kurtosis) and interpret parameters
derive expectation, variance, and the (log-)likelihood
sample from the distribution using NumPy only
visualize PDF/CDF and Monte Carlo simulations
use
scipy.stats.fatiguelifeforpdf,cdf,rvs,fit
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
from scipy import stats, special
from scipy.stats import fatiguelife as fatiguelife_sp
# Plotly rendering (CKC convention)
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
# Reproducibility
rng = np.random.default_rng(7)
np.set_printoptions(precision=6, suppress=True)
SQRT2 = np.sqrt(2.0)
SQRT2PI = np.sqrt(2.0 * np.pi)
1) Title & classification#
Name:
fatiguelifeType: continuous distribution
Standard support: \(x \in (0, \infty)\)
SciPy parameter space (
scipy.stats.fatiguelife(c, loc=0, scale=1)):shape: \(c > 0\) (often written \(\alpha\))
location: \(\mathrm{loc} \in \mathbb{R}\)
scale: \(\mathrm{scale} > 0\) (often written \(\beta\))
Support with
loc/scale: \(x \in (\mathrm{loc}, \infty)\)
We’ll use \((c,\beta)\) for the Birnbaum–Saunders form and note that in SciPy: \(\beta = \mathrm{scale}\).
2) Intuition & motivation#
What it models#
The fatigue-life distribution was introduced to model lifetimes of materials under repeated stress cycles. Intuitively, a material accumulates microscopic damage over many cycles; failure happens when accumulated damage crosses a threshold. This leads to positive support and typically right-skewed lifetimes.
Typical real-world use cases#
Reliability engineering: time-to-failure (or number of cycles to failure) under fatigue
Materials science: crack initiation/propagation times
Survival analysis: positive, skewed durations where lognormal is a competitor
Relations to other distributions#
Normal connection: the transform $\(Z = \frac{1}{c}\left(\sqrt{\frac{X}{\beta}} - \sqrt{\frac{\beta}{X}}\right)\)$ is standard normal.
Lognormal-like behavior: for small \(c\), the distribution concentrates near \(\beta\) and can look approximately lognormal.
Reciprocal symmetry: if \(X\sim\mathrm{BS}(c,\beta)\) then \(\beta^2/X\) has the same distribution. A nice consequence is that the median is exactly \(\beta\).
3) Formal definition#
We write \(X \sim \mathrm{BS}(c,\beta)\) (Birnbaum–Saunders) with \(c>0\) and \(\beta>0\). Let \(\Phi\) and \(\varphi\) be the standard normal CDF and PDF.
CDF#
For \(x>0\):
PDF#
For \(x>0\):
Location–scale form (SciPy)#
SciPy uses the standard rv_continuous transformation:
So for \(x>\mathrm{loc}\), with \(y=(x-\mathrm{loc})/\mathrm{scale}\):
def _validate_params(c: float, scale: float) -> None:
if not (np.isfinite(c) and c > 0):
raise ValueError("c must be finite and > 0")
if not (np.isfinite(scale) and scale > 0):
raise ValueError("scale must be finite and > 0")
def fatiguelife_pdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
'''Fatigue-life PDF in SciPy's loc/scale parameterization.
Notes
-----
Uses the standard form with β=1 on y=(x-loc)/scale and applies the Jacobian 1/scale.
'''
_validate_params(float(c), float(scale))
x = np.asarray(x, dtype=float)
y = (x - loc) / scale
out = np.zeros_like(y, dtype=float)
mask = y > 0
if not np.any(mask):
return out
y_m = y[mask]
s = np.sqrt(y_m)
z = (s - 1.0 / s) / c
phi = np.exp(-0.5 * z**2) / SQRT2PI
jac = (s + 1.0 / s) / (2.0 * c * y_m)
out[mask] = (phi * jac) / scale
return out
def fatiguelife_logpdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
'''Log-PDF (useful for numerical stability).'''
_validate_params(float(c), float(scale))
x = np.asarray(x, dtype=float)
y = (x - loc) / scale
out = np.full_like(y, fill_value=-np.inf, dtype=float)
mask = y > 0
if not np.any(mask):
return out
y_m = y[mask]
s = np.sqrt(y_m)
z = (s - 1.0 / s) / c
log_phi = -0.5 * z**2 - np.log(SQRT2PI)
log_jac = np.log(s + 1.0 / s) - np.log(2.0 * c * y_m)
out[mask] = log_phi + log_jac - np.log(scale)
return out
def fatiguelife_cdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
'''Fatigue-life CDF in SciPy's loc/scale parameterization.'''
_validate_params(float(c), float(scale))
x = np.asarray(x, dtype=float)
y = (x - loc) / scale
out = np.zeros_like(y, dtype=float)
mask = y > 0
if not np.any(mask):
return out
y_m = y[mask]
s = np.sqrt(y_m)
z = (s - 1.0 / s) / c
out[mask] = special.ndtr(z)
return out
def fatiguelife_sf(x, c: float, loc: float = 0.0, scale: float = 1.0):
'''Survival function 1-CDF, computed stably via log_ndtr when possible.'''
_validate_params(float(c), float(scale))
x = np.asarray(x, dtype=float)
y = (x - loc) / scale
out = np.ones_like(y, dtype=float)
mask = y > 0
if not np.any(mask):
return out
y_m = y[mask]
s = np.sqrt(y_m)
z = (s - 1.0 / s) / c
# For large z, 1 - ndtr(z) loses precision; use log_ndtr(-z).
out[mask] = np.exp(special.log_ndtr(-z))
return out
4) Moments & properties#
For \(X\sim\mathrm{BS}(c,\beta)\) with \(c>0\) and \(\beta>0\):
Mean, variance, skewness, kurtosis#
Mean $\(\mathbb E[X] = \beta\left(1 + \frac{c^2}{2}\right).\)$
Variance $\(\mathrm{Var}(X)=\beta^2\,c^2\left(1+\frac{5c^2}{4}\right).\)$
Skewness $\(\gamma_1 = \frac{4c\,(11c^2+6)}{(5c^2+4)^{3/2}}.\)$
Excess kurtosis $\(\gamma_2 = \frac{6c^2\,(93c^2+40)}{(5c^2+4)^2}.\)$
(Scale and location do not change skewness/kurtosis.)
MGF and characteristic function#
A simple closed form is not available. A useful representation uses the normal construction. Let \(Z\sim\mathcal N(0,1)\) and define
Then for SciPy’s parameters \(X=\mathrm{loc}+\mathrm{scale}\,Y\):
MGF (when it exists): $\(M_X(t)=\mathbb E[e^{tX}] = e^{t\,\mathrm{loc}}\,\mathbb E\left[e^{t\,\mathrm{scale}\,Y}\right].\)\( Because \)Y\( grows like \)\approx c^2Z^2\( for large positive \)Z\(, the MGF exists at least for \)\(t < \frac{1}{2\,\mathrm{scale}\,c^2}.\)$
Characteristic function (always exists): $\(\varphi_X(\omega)=\mathbb E[e^{i\omega X}].\)$
Entropy#
The differential entropy \(H(X)=-\mathbb E[\log f(X)]\) has no common elementary closed form.
In practice you compute it numerically (e.g. via scipy.stats.fatiguelife.entropy or Monte Carlo).
def fatiguelife_moments(c: float, loc: float = 0.0, scale: float = 1.0) -> dict:
'''Closed-form mean/var/skew/excess-kurtosis for SciPy's parameters.'''
_validate_params(float(c), float(scale))
mean = loc + scale * (1.0 + 0.5 * c**2)
var = (scale**2) * (c**2) * (1.0 + 1.25 * c**2)
skew = 4.0 * c * (11.0 * c**2 + 6.0) / ((5.0 * c**2 + 4.0) ** 1.5)
excess_kurt = 6.0 * c**2 * (93.0 * c**2 + 40.0) / ((5.0 * c**2 + 4.0) ** 2)
median = loc + scale * 1.0 # median of BS(c,1) is 1
return {
'mean': mean,
'var': var,
'skew': skew,
'excess_kurt': excess_kurt,
'median': median,
}
def normal_expectation_gauss_hermite(func, n: int = 60) -> float:
'''Approximate E[func(Z)] for Z~N(0,1) via Gauss–Hermite quadrature (NumPy-only).'''
xs, ws = np.polynomial.hermite.hermgauss(n)
# For standard normal: E[g(Z)] = 1/sqrt(pi) * sum w_i * g(sqrt(2) * x_i)
zs = SQRT2 * xs
return float(np.sum(ws * func(zs)) / np.sqrt(np.pi))
def fatiguelife_from_normal(z, c: float, loc: float = 0.0, scale: float = 1.0):
'''Deterministic transform mapping N(0,1) samples -> fatigue-life samples.'''
_validate_params(float(c), float(scale))
z = np.asarray(z, dtype=float)
delta = 0.5 * c * z
y = delta + np.sqrt(delta**2 + 1.0)
return loc + scale * (y**2)
def fatiguelife_mgf_gh(t: float, c: float, loc: float = 0.0, scale: float = 1.0, n: int = 80) -> float:
'''MGF via Gauss–Hermite quadrature (best for moderate t).'''
def g(z):
x = fatiguelife_from_normal(z, c=c, loc=0.0, scale=scale)
return np.exp(t * x)
return np.exp(t * loc) * normal_expectation_gauss_hermite(g, n=n)
def fatiguelife_cf_gh(omega: float, c: float, loc: float = 0.0, scale: float = 1.0, n: int = 80) -> complex:
'''Characteristic function via Gauss–Hermite quadrature.'''
def g(z):
x = fatiguelife_from_normal(z, c=c, loc=0.0, scale=scale)
return np.exp(1j * omega * x)
return np.exp(1j * omega * loc) * normal_expectation_gauss_hermite(g, n=n)
params = dict(c=1.0, loc=0.0, scale=2.0)
closed = fatiguelife_moments(**params)
scipy_stats = fatiguelife_sp.stats(params['c'], loc=params['loc'], scale=params['scale'], moments='mvsk')
closed, scipy_stats
({'mean': 3.0,
'var': 9.0,
'skew': 2.5185185185185186,
'excess_kurt': 9.851851851851851,
'median': 2.0},
(3.0, 9.0, 2.5185185185185186, 9.851851851851851))
5) Parameter interpretation#
Shape parameter \(c\) (= c)#
Controls spread, right tail weight, and skewness.
As \(c\to 0\), the distribution collapses toward its median (\(\beta\) in the Birnbaum–Saunders form).
Larger \(c\) produces heavier right tails and larger skewness/kurtosis.
Scale parameter \(\beta\) (= scale)#
Multiplies the variable: if \(Y\sim\mathrm{BS}(c,1)\) then \(X=\beta Y\sim\mathrm{BS}(c,\beta)\).
The median equals \(\beta\).
Location parameter loc#
Shifts the support to \((\mathrm{loc},\infty)\).
In reliability contexts, you sometimes interpret
locas a minimum/threshold lifetime.
x = np.linspace(1e-3, 12, 800)
param_sets = [
(0.2, 0.0, 1.0, 'c=0.2, scale=1'),
(0.5, 0.0, 1.0, 'c=0.5, scale=1'),
(1.0, 0.0, 1.0, 'c=1.0, scale=1'),
(2.0, 0.0, 1.0, 'c=2.0, scale=1'),
]
fig = go.Figure()
for c, loc, scale, label in param_sets:
fig.add_trace(go.Scatter(x=x, y=fatiguelife_pdf(x, c=c, loc=loc, scale=scale), mode='lines', name=label))
fig.update_layout(
title='Fatigue-life PDF: effect of the shape parameter c',
xaxis_title='x',
yaxis_title='density',
width=900,
height=420,
)
fig
6) Derivations#
Expectation#
Start from the constructive representation with \(Z\sim\mathcal N(0,1)\) and \(\delta = \tfrac{cZ}{2}\):
Take expectations. The last term is an odd function of \(Z\) (because \(\delta\) is odd and \(\sqrt{\delta^2+1}\) is even), so it integrates to zero under a symmetric normal.
Finally, for SciPy parameters \(X=\mathrm{loc}+\mathrm{scale}\,Y\):
Variance#
Write \(Y=A+B\) with
\(A = 1+2\delta^2\) (even in \(Z\))
\(B = 2\delta\sqrt{\delta^2+1}\) (odd in \(Z\))
Then
and \(\mathbb E[AB]=0\) because \(AB\) is odd in \(Z\). Using \(\mathbb E[Z^2]=1\) and \(\mathbb E[Z^4]=3\), you get
So
Scaling gives
Likelihood (i.i.d. sample)#
Given data \(x_1,\dots,x_n\) and parameters \((c,\mathrm{loc},\mathrm{scale})\), define \(y_i=(x_i-\mathrm{loc})/\mathrm{scale}\). Validity requires \(c>0\), \(\mathrm{scale}>0\), and all \(y_i>0\).
For each \(i\):
where \(z_i=\tfrac{1}{c}(\sqrt{y_i} - 1/\sqrt{y_i})\). The total log-likelihood is the sum over \(i\).
def fatiguelife_nll(params, data):
'''Negative log-likelihood for (c, loc, scale).'''
c, loc, scale = params
if not (np.isfinite(c) and c > 0 and np.isfinite(scale) and scale > 0 and np.isfinite(loc)):
return np.inf
y = (data - loc) / scale
if np.any(y <= 0):
return np.inf
return -float(np.sum(fatiguelife_logpdf(data, c=c, loc=loc, scale=scale)))
# Quick check: compare our logpdf to SciPy's for random x>0
c0, loc0, scale0 = 0.8, 0.0, 1.7
x_test = rng.uniform(1e-4, 8, size=10_000)
err = np.max(
np.abs(
fatiguelife_logpdf(x_test, c=c0, loc=loc0, scale=scale0)
- fatiguelife_sp.logpdf(x_test, c0, loc=loc0, scale=scale0)
)
)
err
1.1368683772161603e-12
7) Sampling & simulation (NumPy only)#
The normal-construction gives a very clean sampler.
Algorithm#
To sample \(X\sim\mathrm{BS}(c,\beta)\) (with optional loc):
Sample \(Z\sim\mathcal N(0,1)\).
Compute \(\delta = \tfrac{cZ}{2}\).
Compute \(Y = \delta + \sqrt{\delta^2+1}\).
Return \(X = \mathrm{loc} + \beta\,Y^2\).
This is vectorized and uses only numpy.random + elementary operations.
def fatiguelife_rvs_numpy(c: float, loc: float = 0.0, scale: float = 1.0, size=1, rng=None):
'''Sample from the fatigue-life distribution using NumPy only.'''
_validate_params(float(c), float(scale))
if rng is None:
rng = np.random.default_rng()
z = rng.standard_normal(size=size)
return fatiguelife_from_normal(z, c=c, loc=loc, scale=scale)
samples_np = fatiguelife_rvs_numpy(1.0, loc=0.0, scale=1.5, size=5, rng=rng)
samples_np
array([2.243943, 1.852748, 0.772983, 7.725928, 0.481773])
8) Visualization#
We’ll compare:
the theoretical PDF/CDF
Monte Carlo samples (histogram + empirical CDF)
c0, loc0, scale0 = 1.0, 0.0, 2.0
n = 80_000
samps = fatiguelife_rvs_numpy(c0, loc=loc0, scale=scale0, size=n, rng=rng)
x_grid = np.linspace(1e-3, np.quantile(samps, 0.995), 600)
# Empirical CDF
xs_sorted = np.sort(samps)
ecdf = np.arange(1, n + 1) / n
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=('PDF + Monte Carlo histogram', 'CDF + empirical CDF'),
)
# PDF panel
hist_y, hist_x = np.histogram(samps, bins=120, density=True)
hist_xc = 0.5 * (hist_x[1:] + hist_x[:-1])
fig.add_trace(go.Bar(x=hist_xc, y=hist_y, name='MC histogram', opacity=0.45), row=1, col=1)
fig.add_trace(
go.Scatter(
x=x_grid,
y=fatiguelife_pdf(x_grid, c=c0, loc=loc0, scale=scale0),
mode='lines',
name='theoretical pdf',
line=dict(width=3),
),
row=1,
col=1,
)
# CDF panel
fig.add_trace(
go.Scatter(
x=x_grid,
y=fatiguelife_cdf(x_grid, c=c0, loc=loc0, scale=scale0),
mode='lines',
name='theoretical cdf',
line=dict(width=3),
),
row=1,
col=2,
)
fig.add_trace(
go.Scatter(
x=xs_sorted[::40],
y=ecdf[::40],
mode='markers',
name='empirical cdf',
marker=dict(size=4),
),
row=1,
col=2,
)
fig.update_layout(width=1000, height=420, legend=dict(orientation='h', y=1.15))
fig.update_xaxes(title_text='x', row=1, col=1)
fig.update_yaxes(title_text='density', row=1, col=1)
fig.update_xaxes(title_text='x', row=1, col=2)
fig.update_yaxes(title_text='F(x)', row=1, col=2)
fig
9) SciPy integration (scipy.stats.fatiguelife)#
SciPy’s object supports the standard continuous-distribution API:
pdf,logpdf,cdf,sfrvsfor samplingfitfor parameter estimationstats,moment,entropy, …
The call signature is:
scipy.stats.fatiguelife(c, loc=0, scale=1)
where c is the shape parameter and scale corresponds to \(\beta\).
# Sampling via SciPy
samps_sp = fatiguelife_sp.rvs(c0, loc=loc0, scale=scale0, size=n, random_state=rng)
# Fit parameters (returns: c_hat, loc_hat, scale_hat)
fit_c, fit_loc, fit_scale = fatiguelife_sp.fit(samps_sp)
fit_c, fit_loc, fit_scale
(0.9997818369711162, 0.001552102122452088, 2.0038068173424053)
# Compare fitted PDF to the true PDF
x_grid = np.linspace(1e-3, np.quantile(samps_sp, 0.995), 600)
fig = go.Figure()
fig.add_trace(go.Histogram(x=samps_sp, nbinsx=120, histnorm='probability density', name='data', opacity=0.45))
fig.add_trace(
go.Scatter(
x=x_grid,
y=fatiguelife_sp.pdf(x_grid, c0, loc=loc0, scale=scale0),
mode='lines',
name='true pdf',
line=dict(width=3),
)
)
fig.add_trace(
go.Scatter(
x=x_grid,
y=fatiguelife_sp.pdf(x_grid, fit_c, loc=fit_loc, scale=fit_scale),
mode='lines',
name='fitted pdf',
line=dict(width=3, dash='dash'),
)
)
fig.update_layout(
title='SciPy fit: true vs fitted PDF',
xaxis_title='x',
yaxis_title='density',
width=900,
height=420,
)
fig
10) Statistical use cases#
Hypothesis testing (goodness-of-fit)#
A common approach is the probability integral transform (PIT). If the model is correct and parameters are known, then
So we can test uniformity of \(U_i\) using a KS test.
Important caveat: if you fit parameters on the same data and then run KS, the p-values are not exact. A more principled approach is a parametric bootstrap.
Bayesian modeling#
For positive parameters (\(c>0\), \(\beta>0\)) it’s common to model them on a log-scale, e.g.
Then the posterior is proportional to
Below we’ll do a simple grid posterior demonstration (no external Bayesian libraries).
Generative modeling#
fatiguelife is a convenient generator for positive, skewed data.
For example, it can model synthetic failure times in simulations, or serve as a base distribution in mixture/latent-variable models.
# Goodness-of-fit via PIT + KS (illustrative)
# We'll use the fitted SciPy params from above.
U = fatiguelife_sp.cdf(samps_sp, fit_c, loc=fit_loc, scale=fit_scale)
ks = stats.kstest(U, 'uniform')
ks
KstestResult(statistic=0.0016533666280400539, pvalue=0.9807346694212495, statistic_location=0.7203466333719599, statistic_sign=1)
# Simple Bayesian grid posterior for (c, scale) with loc fixed at 0
# (This is a pedagogical demo; for serious work use MCMC/VI.)
loc_fixed = 0.0
# Keep a smaller dataset for speed
data = samps_sp[:2000]
# Grid over log-parameters
logc_grid = np.linspace(np.log(0.2), np.log(3.0), 120)
logs_grid = np.linspace(np.log(0.3), np.log(6.0), 120)
# Log-likelihood on the grid (vectorized over scale for each c)
ll = np.zeros((logc_grid.size, logs_grid.size), dtype=float)
scale_grid = np.exp(logs_grid)
for i, c in enumerate(np.exp(logc_grid)):
ll[i, :] = np.sum(fatiguelife_sp.logpdf(data[:, None], c, loc=loc_fixed, scale=scale_grid[None, :]), axis=0)
# Priors: log c ~ N(0, 1^2), log scale ~ N(0, 1^2)
log_prior = -0.5 * (logc_grid[:, None] ** 2 + logs_grid[None, :] ** 2)
log_post_unnorm = ll + log_prior
log_post = log_post_unnorm - np.max(log_post_unnorm)
post = np.exp(log_post)
post /= np.sum(post)
# Posterior summaries
c_mean = float(np.sum(np.exp(logc_grid)[:, None] * post))
scale_mean = float(np.sum(np.exp(logs_grid)[None, :] * post))
(c_mean, scale_mean)
(1.0045700659359345, 1.9890536073446352)
# Visualize the posterior over (c, scale)
fig = go.Figure(
data=go.Heatmap(
x=np.exp(logs_grid),
y=np.exp(logc_grid),
z=post,
colorscale='Viridis',
colorbar=dict(title='posterior mass'),
)
)
fig.update_layout(
title='Grid posterior p(c, scale | data) with log-normal priors (loc fixed at 0)',
xaxis_title='scale (β)',
yaxis_title='c',
width=850,
height=520,
)
# Add markers for truth and posterior mean
fig.add_trace(go.Scatter(x=[scale0], y=[c0], mode='markers', name='true', marker=dict(size=10, symbol='x')))
fig.add_trace(
go.Scatter(
x=[scale_mean],
y=[c_mean],
mode='markers',
name='posterior mean',
marker=dict(size=10, symbol='circle-open'),
)
)
fig
11) Pitfalls#
Invalid parameters: must have
c > 0,scale > 0. Withloc, all observations must satisfy \(x>\mathrm{loc}\).Parameterization confusion: many texts use \((\alpha,\beta)\) for (shape, scale); SciPy uses
c(shape) andscale(\(=\beta\)).Numerical stability:
Use
logpdfinstead ofpdffor likelihood work.For extreme right tails, prefer
sf/logsfinstead of1-cdf.
MGF domain: the MGF exists only for sufficiently small positive \(t\); avoid relying on it outside its radius of convergence.
Goodness-of-fit after fitting: KS p-values are not exact when parameters are estimated from the same sample; consider bootstrap.
12) Summary#
fatiguelife(Birnbaum–Saunders) is a continuous distribution on \((0,\infty)\) used for fatigue failure times.It has a clean normal connection via the transform \(\tfrac{1}{c}(\sqrt{x/\beta}-\sqrt{\beta/x})\).
Mean/variance/skewness/kurtosis are available in closed form; the median equals \(\beta\).
Sampling is easy with a NumPy-only transformation from a standard normal.
In practice,
scipy.stats.fatiguelifeprovidespdf,cdf,rvs,fit, and numericalentropy.